Homotopy techniques for multiplication 
modulo triangular sets 



Alin Bostan (alin.bostan@inria.fr) 



Algorithms Project, INRIA Rocquencourt, 78153 Le Chesnay Cedex, France 



Muhammad Chowdhury (mchowdh3@csd.uwo.ca) 

Computer Science Department, The University of Western Ontario, London, Ontario, Canada 

Joris van der Hoeven (Joris.Vanderhoeven@math.u-psud.fr) 

CNRS, Departement de Mathematiques, Universite Paris-Sud, 91405 Orsay Cedex, France 



Computer Science Department, The University of Western Ontario, London, Ontario, Canada 
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families of examples, for which no such result was known before. Applications are given to 
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1. Introduction 

Triangular families of polynomials are a versatile data structure, well adapted to en- 
code geometric problems with some form of symmetry [3J [TSJ E3J EH] ■ However, in spite 
of this, many complexity questions are still not answered in a satisfying manner. 
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A high-level question is to provide sharp estimates on the cost of solving polynomial 
systems by means of triangular representations. This problem has a geometric nature; it 
itself relies on several difficult lower-level questions, such as the cost of basic operations 
with triangular sets. In this paper, we address one such question: the arithmetic cost of 
multiplication of polynomials modulo a triangular set. This justifiably stands as a central 
question, since many higher-level routines are built on top of it, such as inversion [2011551 
[22] . lifting techniques [TU] for modular algorithms, solving systems of equations [2T| . etc. 



Problem statement, overview of our results. We work in R[Xi, 
and we are given a set of relations of the form 



T n (X\, . . . , X n ) 



, X n ], where R is a ring, 



T 2 (X U X 2 ) 

The polynomials T form a triangular set: for all i, Ti is in R[Xi, . . . , XJ, is monic in 
Xi and reduced modulo (Tx, . . . , in the sense that deg(Ti,Xj) < deg(Tj, Xj) for 

j < i. As an aside, note that the case where Ti is not monic but with a leading coefficient 
invertible modulo (Xi, . . . , Xj_i) reduces in principle to the monic case; however, inversion 
modulo (Xi, . . . , Ti_i) remains a difficult question [llj . of complexity higher than that of 
multiplication. 

As input, we consider two polynomials A, B reduced modulo (T). The direct approach 
to multiply them modulo (T) is to perform a polynomial multiplication, followed by the 
reduction modulo (T), by a generalization of Euclidean division. As far as complexity 
is concerned, when the number of variables grows, this kind of approach cannot give 
linear time algorithms. Consider for instance the case where all Ti have degree 2 in their 
main variables Xi. Then, A and B both have 2™ monomials, but their product before 
reduction has 3 n monomials; after reduction, the number of monomials is 2™ again. If we 
let 5 = 2™ be a measure of the input and output size, the cost of such an algorithm is at 
least 3" = <5 log 2 3 . 

In this paper, we show that a different approach can lead to a quasi-linear time algo- 
rithm, in cases where the monomial support of T is sparse, or when the polynomials in 
T have a low total degree. This will for example be the case for systems of the form 



X n ~ 2X n _i 



X<2 — *2<X\ 

xl 



or 



Xn ~ X n -1 



xl- 
xl 



X, 



(1) 



whose applications are described later on. Our result also applies to the following con- 
struction: start from F E R[X], say F — X 3 — X 2 + X — 3, and define the so-called 
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'Cauchy modules" [27], which are in effective Galois theory [271 HI [25]: 



X 3 + X 2 + X l -l 

Xf + X 2 X l -X 2 + X 2 1 -X l + l (2) 
Xl - X\ + Xi - 3. 

For examples (1) and (2), our algorithms give the following results: 

• for T as in (1), multiplication modulo (T) can be performed in quasi-linear time 0~(5), 
where S = 2™ is the input and output size, and where 0~(S) stands for 0(S(logS) ^ 1 '). 

• for T as in (2), with n — deg(F), multiplication modulo (T) can be performed in 
quasi-linear time 0~(S), where 5 = n\ is the input and output size. 

No previous algorithm was known featuring such complexity estimates. 

Our approach. To obtain this quasi-linear cost, we have to avoid multiplying A and B 
as polynomials. Our solution is to use evaluation and interpolation techniques, just as 
FFT multiplication of univariate polynomials is multiplication modulo X n — 1 . 

Fast evaluation and interpolation may not be possible directly, if T does not have roots 
in R (as in the previous examples). However, they become possible using deformation 
techniques: we construct a new triangular set U with all roots in R, and multiply A and 
B modulo S = 77T + (1 — J7)U, where r\ is a new variable. The triangular set S has roots in 
R[W]> by Hensel's lemma, so one can use evaluation-interpolation techniques over R[[r]]]. 

This idea was introduced in [22] , but was limited to the case where all polynomials in 
T are univariate: Tj was restricted to depend on Xj only, so this did not apply to the 
examples above. Here, we extend this idea to cover such examples; our main technical 
contribution is a study of precision-related issues involved in the power series computa- 
tions, and how they relate to the monomial support of T. 

Previous work. It is only recently that fast algorithms for triangular representations have 
been throughly investigated; thus, previous results on efficient multiplication algorithms 
are scarce. All natural approaches introduce in their cost estimate an overhead of the 
form k n , for some constant k. 

The main challenge (still open) is to get rid of this exponential factor unconditionnally: 
we want algorithms of cost 0~(5), where 5 = d\ ■ ■ ■ d n is the number of monomials in A, 
B and their product modulo (T). For instance, with = Xf*, the first complexity result 
of the form 0(<5 1+e ), for any e > 0, was in [55] . 

The previous work [22] gives a general algorithm of cost 0~(4 n 5). That algorithm 
uses fast Euclidean division; for polynomials of low degree (e.g., for di = 2), the naive 
approach for Euclidean division can actually give better results. Previous mentions of 
such complexity estimates (with a constant higher than 4) are in [20] . 

As said above, in [22], one also finds the precursor of the algorithm presented here; 
the algorithm of [22] applies to families of polynomials having Tj € R[Xj], and achieves 
the cost 0(5 1+e ) for any e > 0. In that case, the analysis of the precision in power series 
computation was immediate. Our main contribution here is to perform this study in the 
general case, and to show that we can still achieve similar costs for much larger families 
of examples. 



^3(Xi, X2, X3) 

F 2 (Xi, X 2 ) 
FAXi) 



F 2 (X U X 2 )~F(X 1 ,X 3 ) 
X2 — X3 

F 1 (X 1 )-F 1 (X 2 ) 

Xl —X-2 

F(Xi) 
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Basic notation. Let d = {d\, ■ ■ ■ ,d n ) be a vector of positive integers. In what follows, 
these will represent the main degrees of the polynomials in our triangular sets; without 
loss of generality, we will thus always suppose di > 2 for all i. 

Recall that R is our base ring and that X\, . . . ,X n are interminates over R. We let 
M d be the set of monomials 

M d = {XI 1 -~X%> \ Q<e l <d i for alH }. 

We denote by Span(Ma) the free R-submodule generated by Md in R[-Xi, ■ • ■ ,X n ]: 

Span(Md) = { A = a m m \ a m £ R for all m £ Md }. 

mEMd 

This is thus the set of polynomials A in R[Xi, . . . , X n ], such that deg(Ai, Xi) < di holds 
for all i. Finally, we let <5d be the product 5 d = d± ■ ■ ■ d n ; this is the cardinality of Md. 
Remark that since all di are at least 2, we have the bounds 

2 n <S d and ^ f d 1 ---d i < 25 d . 

The former plainly follows from the inequality 2 < d^ the latter comes from observing 
that d\ ■ ■ ■ di2 n ~ l < d\- ■ ■ d n = S d ', this yields d\ ■ ■ ■ di < 5 d /2 n ~' 1 , from which the claim 
follows by summation. 

The multi-degree of a triangular set T = (Ti, . . . , T n ) is the n-uple d = (d\, . . . , d n ), 
with di = deg(T i , XA 1<i<n . In this case, R[Xi, . . . , X n ]/ (T) is a free R-module isomorphic 
to Span(Afd)- We say that a polynomial A € R[-X"i, ■ ■ ■ ,X n ] is reduced with respect to 
T if deg(A,Xi) < di holds for all i. For any A s R[-Yi, ■ ■ • ,X n ), there exists a unique 
A' 6 R[-X"i, . . . , X n ], reduced with respect to T, and such that A — A' in is the ideal (T). 
We call it the normal form of A and write A' = A mod (T) . 

Outlook of the paper. In Section 2, we introduce some basic complexity notation. The 
next section presents basic evaluation-interpolation algorithms for so-called equipro- 
jectable sets, which arc extensions of algorithms known for univariate polynomials. We 
deduce our multiplication algorithm in Section 4; examples, applications and experimen- 
tal results are in Sections 5 and G. 

2. Preliminaries 

Big-0 notation is delicate to use in our situation, since our estimates may depend on 
several (possibly an unbounded number of) parameters (typically, the multi-degree of our 
triangular sets). Hence, whenever we use a big-0 inequality such as f € O(g), it is implied 
that there exists a universal constant A such that f(vi, . . . ,v s ) < Xg(v\, . . . ,v s ) holds for 
all possible values of the arguments. When needed, we use explicit inequalities. Finally, 
the notation / € 0~(g) means that there exists a constant a such that / £ 0(glog(g) a ), 
where the big-0 is to be understood as above. 

Our complexity estimates count additions, multiplications, and inversions, when they 
are possible. We denote by M : N — > N a function such that over any ring, polynomials 
of degree less than d can be multiplied in M(d) operations, and which satisfies the super- 
linearity conditions of [TH Chapter 8]. Using the algorithm of Cantor-Kaltofen [S], one 
can take M(d) £ 0(dlg(d) lglg(ef)), with lg(d) = log 2 max(d, 2). 

We next let Co : N — > N be a function such that Co(d) > d holds for all d and such 
that we have, over any ring R: 
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(1) for any xi, . . . , Xd in R and any polynomial A E R[X] of degree less than d, one can 
compute all values A(xi) in Co(d) additions and multiplications in R; 

(2) for any x\, . . . ,Xd in R, one can compute the coefficients of the polynomial (X — 
Xi) ■ ■ ■ (X — x<i) in Co(<i) additions and multiplications in R; 

(3) for any x±, . . . , x& in R, with x^ — Xj a unit for i ^ j, and any values vi, . . . , in R, 
one can compute the unique polynomial A £ R[X] of degree less than d such that 
A(x{) = Vi holds for all i in Co(c?) operations in R. 

By the results of jT5j Chapter 10], one can take Co(d) 6 0(M(d) log(d)). We continue 
with the well-known fact that the function Co also enables us to estimate the cost of 
lifting power series roots of a bivariate polynomial by Newton iteration. In the following 
lemma, 77 is a new variable over R. 

Lemma 1 For any polynomial T in R[rj, X], monic in X and with deg(T, X) = d, if the 
roots o/T(0,X) are known and have multiplicity 1, one can compute the roots of T{ri, X) 
in R[?7]/(?/) in 0(Co(rf)M(£)) operations in R. 

Proof. The algorithm consists in lifting all roots of T in parallel using Newton iteration, 
using fast evaluation to compute the needed values of T and dT/dX; it is given in 
Figure 1, where we use a subroutine called EvalUnivariate to do the evaluation. Each 
pass through the loop at line 2 takes two evaluations in degree d and d inversions, with 
coefficients that are power series of precision £'. Using [15] Chapter 9], the cost is thus 
2Co(d)M(£') + XdM(£'), for some constant A. Using the super-linearity of the function M, 
the conclusion follows. □ 

Remark 1 When performing all multiplications in R[?7, X\j (rf) using Kronecker's method, 
a more precise cost analysis yields the bound 0(M(d£) log(<i)) instead of 

O(C (d)M(£)) = 0(M(d)M(£)logd). 

If £ — O(d), then we usually have M(d£) = 0(M(d)£), which makes the new bound 
slightly better. However, this improvement only has a minor impact on what follows, so 
it will be more convenient to use the technically simpler bound from the lemma. 



LiftRoots(T,ai, . . . ,a d ,£) 

1 £' <- 2 

2 while f < £ do 

2.1 vi, . . . , Vd *— EvalUnivariate(T, 01, . . . , ad) mod tj 1 

2.2 wi, ...,Wd*— EvalUnivariate(9T/9X, ai, . . . , ad) mod rj 

2.3 for i = 1, . . . , d do 

2.3.1 ai <— ai — Vi/wi mod rf 

2.4 £' <- 21' 

3 return [ai mod X £ \ 1 < i < d] 



Fig. 1. Lifting all roots of a univariate polynomial. 

To obtain simpler estimates, we let C(d) = ACo(d), where A > 1 is the constant implied 
in the big-0 estimate in the former lemma. Hence, problems (1), (2) and (3) above can 
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be dealt with in C(d) operations, and the lifting problem of the previous lemma can be 
solved in C(d)M(£) operations. Finally, we introduce another short-hand notation: for a 
multi-degree d = (d\, . . . , d n ), we write 

L(d)=E^<n^ (3) 

with d = maxi<„G?i. In view of the estimates on C, we also have the upper bound 
L(d) € 0(lg(6d) 3 ), which shows that L(d) is of polylogarithmic growth in Sd- 



3. Evaluation and interpolation at equiprojectable sets 

In this section, we recall from [3] the definition of equiprojectable sets. We prove that 
one can perform evaluation and interpolation at, and construct the vanishing ideal of, 
equiprojectable sets in linear time, up to logarithmic factors. We deduce an algorithm for 
multiplication modulo the vanishing ideal of such sets with a similar complexity. These 
results extend those given in pHl 122) . which dealt with the case of points on a regular 
grid. The extension to our more general context is rather straightforward, but to our 
knowledge, it has not appeared in print before. 

In all this section, R is a ring; we study subsets of R n and their successive projections 
on the subspaces R l , for i < n. For defmiteness, we let R° be a one-point set. Then, for 
1 < j < i < we let %i d be the projection 

n id : R' -> R' 

(an,..., or*) h-> (xx,...,Xj); 

if j = 0, we adapt this definition by letting 71^0 be the constant map R l — > R°. Finally, 

since this is the projection we use most, we simply write 7r = ~K n>n —i for the projection 
R« _> R n-1 

If V is a subset of R", for (3 in 7r(V), we let Vp be the fiber V n tt" 1 ^). Hence, if (3 
has coordinates (fix, ■ ■ ■ , /3 n _i), the points in Vp have the form (/3i, . . . , a), for some 
values a in R. In all that follows, a finite set is by convention non-empty. 

3.1. Equiprojectable sets 

Let V be a finite set in R™. Equiprojectability is a property of V that describes a 
combinatorial regularity in the successive projections of V. For n = 0, we say that 
the unique non-empty subset of R° is equiprojectable. Then, for n > 0, V C R™ is 
equiprojectable if the following holds: 

• the projection n(V) is equiprojectable in R" _1 , and 

• there exists an integer d n such that for all (3 in ir(V), the fiber Vp has cardinality d n . 

The vector d = (d±, . . . ,d n ) is called the multi-degree of V. Remark that for n = 1, 
any finite V C R is equiprojectable. One easily sees that if V is equiprojectable, its 
cardinality equals Sd = di ■ ■ ■ d n ; more generally, Tr 7li i(V) C R 1 is equiprojectable of 
cardinality d\ ■ ■ ■ di. When R is a perfect field, it is proved in [3] that equiprojectable sets 
are exactly the zero-sets of triangular sets that generate radical ideals in R[Xi, . . . ,X n ]; 
we will discuss this in more detail in Subsection 3.4. 
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We give first a slightly more precise notation for the fibers Vp: if V is equiprojectable, 
then for all = (0\, . . . , n -i) G tt(V), there exist exactly d n pairwise distinct values 
vp = [ap t i, . . . , a/3,d„] in R such that 

Vp = . .,/3„_i, a/3,i) | ap }i G vp] 

and thus 

V = [(J3i,..., f3 n -i,a ,i) | = (0i,... ,0 n -i) e 7r(V), e^, 1 < i< d„]. 

For instance, n-dimensional grids are special cases of equiprojectable sets, where vp is 
independent of 0. Remark also that for some special choices of vp, improvements in the 
algorithms below are possible (e.g., for vp translates of points in geometric progression, 
using the algorithm of [2]). 

3.2. Evaluation 

Let V be an equiprojectable set of multi-degree d = (di, . . . ,d n ), and let Aid, <5d be 
as in Section 1. We denote by Evaly the evaluation map 

Evaly : Span(M d ) ->■ R 5d 

F i ► [F(a] | a G V]. 

We let CEvai be a function such that for any V equiprojectable of multidegree d, the map 
Evaly can be evaluated in CEvai(d) operations. In one variable, with n = 1 and d = (di), 
CEvai(d) simply describes the cost of evaluating a polynomial of degree less than di at 
di points of R, so we can take CEvai(d) = C(rfi). More generally, we have the following 
quasi-linear time estimate. 

Proposition 1 One can take CE va i(d) < <5d£(d). 



Proof. We will use a straightforward recursion over the variables X n , . . . , X\ . Let W = 
7r(V), let e = (di, ■ ■ ■ , d n -i) be the multi-degree of W and let A(Xi, . . . , X n ) G Span(Md) 
be the polynomial to evaluate. We write 

A = ^ Ai(Xi, . . . , X n _i)X l n , 

i<d n 

with Ai in Span(A^ e ), and, for in R™ -1 , we define 

Ap= Y,M0)K e R[x n }. 

i<d n 

Hence, for — (0i, . . . , n -i) in R" _1 and x in R, A(0i, . . . ,/3 n _i, x) — Ap(x). As a 
consequence, to evaluate A at V, we start by evaluating all Ai at all points G W. This 
gives all polynomials Ap, which we evaluate at the fibers Vp. The algorithm is given in 
Figure 2. From this, we deduce that we can take CEvai satisfying the recurrence 

CEvai(^i, ■ ■ ■ i d n ) < CEvai(di, ■ ■ • , d n -i) d n + di ■ ■ ■ d n -iC(d n ). 

This implies 

CEva\(di,...,d n ) < Sd ( ~4~ L ~ i 
which proves the proposition. □ 
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Eva\(A,V) 

1 if n = return [A] 

2 W^tt(V) 

3 for i = 0, . . . , d n — 1 do 

3.1 Aj <- coeff(A,X„,i) 

3.2 val[i] <- Eva\(Ai,W) 

/* val[i] has the form [A t ((3) \ (3 G W] */ 

4 for /3 in VF do 

5 return [EvalUnivariate(^4a, u^) | (3 £ W] 



Fig. 2. Evaluation algorithm. 

5.5. Interpolation 

Using the same notation as above, the inverse of the evaluation map is interpolation 
at V: 

lnterp v : R 5d -> Span(M d ) 

[F(a) | a e V] ^ F. 

For this map to be well-defined, we impose a natural condition on the points of V. Let 
W = 7r(V) G R™ _1 . We say that V supports interpolation if 

• if n > 1, W supports interpolation, and 

• for all (3 in W and all x, x' in vp, x — x' is a unit; 

if the base ring is a field, this condition is vacuous. We will see in the following proposition 
that if V supports interpolation, then the map Interpy is well-defined. Moreover, we let 
Qnterp be such that, for V cquiprojectable of multi-degree d, if V supports interpolation, 
then the map lnterp v can be evaluated in C| nte rp(d) operations (including inversions). 
Proposition 2 If V supports interpolation, the map lnterp y is well-defined. Besides, 
one can take Q n terp(d) < 5dL(d). 

Proof. If n = 0, we do nothing; otherwise, we let W = tt(V). The set of values to 
interpolate at V has the shape [f a \ a £ V] e R 5d ; we can thus rewrite it as [fp \ (3 £ W] , 
where each fp is in R d ". 

Since V supports interpolation, for (3 in W, there exists a unique polynomial Ap e 
R[X n ] of degree less than d n , such that Eva\(Ap, vp) = fp. Applying the algorithm re- 
cursively on the coefficients of the polynomials Ap, we can find a polynomial A such 
that A{(3, X n ) = Ap(X n ) holds for all [3 € W. Then, the polynomial A satisfies our con- 
straints. This provides a right-inverse, and thus a two-sided inverse for the map Eva I. The 
algorithm is given in Figure 3; we use a subroutine called I nterp Univariate for univariate 
interpolation. As for evaluation, we deduce that we can take C| nterp satisfying 

C|nterp(^l> • • • i ^n) ^ C| n terp 

(di, d n -i) d n + di ■ ■ ■ d„_iC(d n ), 
which gives our claim, as in the case of evaluation. □ 
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Interp(/,V) 


1 


if n = return [/] 


2 




3 


for /3 in 1U do 


3.1 


<— lnterpUnivariate(//3, vp) 


4 


for i = 0, . . . , d n — 1 do 


4.1 


c, <- [coeff( J 4 i g,X„,i) liSeW] 


4.2 


A 4 <- lnterp(c i7 TU) 


5 


return £ 4<(in j^X* 



Fig. 3. Interpolation algorithm. 



3.4- Associated triangular set 

Next, we associate to an equiprojectable set V C R™ of multi-degree d = (di, . . . ,d n ) 
a triangular set T = (Ti, . . . , T„) of the same multi-degree, which vanishes on V. As soon 
as V supports interpolation, the existence of T is guaranteed (and is established in the 
proof of the next proposition). Uniqueness holds as well: if (Ti, . . . , T n ) and (T{, . . . , T^) 
both vanish on V and have multi-degree d, then for all i, Ti — T[ vanishes at V as well 
and is in Span(Md); hence, it is zero. We call T the associated triangular set; if R is a 
field, T is a lexicographic Grobner basis of the vanishing ideal of V. 
Proposition 3 Given an equiprojectable set V of multi-degree d that supports interpo- 
lation, one can construct the associated triangular set T in time 0(<$ d L(d)). 



Proof. We proceed inductively, and suppose that we already have computed T±, . . . , T„_i 
as the associated triangular set of W = n(V). We will write d = (di,...,d„) and 
e = (di, . . . ,dn_i). 

For /? in W, let Tp be the polynomial H aev0 {X n - a) e R[X n ]. For j < d n , let 
further Tj, n be the polynomial in Span(M e ) that interpolates the jth coefficient of the 
polynomials Tp at W; for j = d n , we take T dn ^ n = 1. We then write T„ = J2j<d n ^j,nX^- 
this polynomial is in R[Xl, . . . , X n ], monic of degree d n in X n , has degree less than 
di in Xi, for i < n, and vanishes on V. Thus, the polynomials T = (Ti, . . . , T„) form 
the triangular set we are looking for. The algorithm is in Figure 4; we use a function 
Poly From Roots to compute the polynomials Tp. 

For a given (3 in W, the function PolyFromRoots computes Tp in C(d„) base ring opera- 
tions; this implies that given T\, . . . , T n _i, one can construct T„ using di • • • d n -iC(d n ) + 
Cinterp(di, ■ ■ ■ , d n _i)d„ operations. The total cost for constructing all Ti is thus at most 

^ dl • • • di-lC(dj) + C|nterp(dl, ■ • • , d;_l)dj. 

Using the trivial bound di • • • dj < <5 d for the left-hand term, and the bound given in 
Proposition 2 for the right-hand one, we get the upper bounds 

z<n i<n J'Si't— 1 z<n z<n 

Using the upper bound J2i<n d\ - • • di < 2<5d, we finally obtain the estimate 3<5dL(d). □ 
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AssociatedTriangularSet(V, n) 

1 if n = return [] 

2 W <- 7r(V) 

3 (Ti, . . . , T n _i) «— AssociatedTriangularSet(W, n) 

4 for in W do 

4.1 T/3 <- Poly From Roots(u )3 ) 

5 for j = 0, . . . , d n — 1 do 

5.1 T ijW «- lnterp([coeff(T^,X n , j) | (3 e W], 

6 return £ Kd „ T s> ^ + 



Fig. 4. Associated triangular set of V . 

3.5. Multiplication 

Using our evaluation and interpolation algorithms, it becomes immediate to perform 
multiplication modulo a triangular set T associated to an equiprojectable set. 
Proposition 4 Let V C R n be an equiprojectable set of multi- degree d = (di, . . . ,d n ) 
that supports interpolation, and let T be the associated triangular set. Then one can 
perform multiplication modulo (T) in time 0(<5d£(d)). 

Proof. The algorithm is the same as in j2SJ Section 2.2], except that we now use the 
more general evaluation and interpolation algorithms presented here. Let A and B be 
reduced modulo (T), and let C = AB mod (T). Then for all a in V, C(a) = A(a)B(a). 
Since C is reduced modulo (T), it suffices to interpolate the values A(a)B(a) to obtain 
C. The cost is thus that of two evaluations, one interpolation, and of all pairwise pairwise 
products; the bounds of Propositions 1 and 2 conclude the proof. □ 



M\A{A,B,V) 

1 VaU <- Eval(A V) 

2 Val B <- Eval(S,y) 

3 Val c «- [ VaU(a)Val B (a) | a G V ] 

4 return lnterp(Valc, V) 



Fig. 5. Multiplication algorithm. 



4. Homotopy techniques for multiplication 

Let T be a triangular set in R[Xi, . . . , X n ], We saw in the previous section that if T 
has all its roots in R, and if V(T) supports interpolation, then multiplication modulo (T) 
can be done in quasi-linear time. In this section, we extend this approach to an arbitrary 
T by setting up an homotopy between T and a new, more convenient, triangular set U. 
This extends the approach of [jJS] Section 2.2], which dealt with the case where is in 
R[Xi] for alii 
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Let d be the multi-degree of T and assume that there exists an equiprojectable set V 
in R™ which supports interpolation and has multi-degree d. Let U be the triangular set 
associated to V and let 77 be a new variable. We then define the set S in R[[r?]] [Xi , . . . , X n ] 
by 

S, = riT, + (1 - n)U it 1 < i < n. 
Since U and T have the same multi-degree d, this set S is triangular, with multi-degree d. 

In Subsection 4.1, we prove that S has all its roots in R[[r]]]. Thus, we can use 
evaluation-interpolation techniques to do multiplication modulo (S); this will in turn 
be used to perform multiplication modulo (T). 

The algorithm involves computing with power series; the quantity that will determine 
the cost of the algorithm will be the required precision in 77. For e = (ei, . . . , e n ) in N", 
we define 

flo(ei, . . . , e„) = deg(X 1 ei ■ • • X^ mod (S),fj) 

and 

H(e 1 ,...,e n )= max H (e' l7 . . . ,e' n ). 

e^<ei, e' n <e n 

Let us then define r = H(2dx — 2, ... , 2d n — 2). Subsection 4.2 shows that multiplication 
modulo (T) can be performed in time 0~{5^r). Finally, in Subsection 4.3, we give upper 
bounds on r that are determined by the monomial support of S; this is the technical core 
of this article. 

4-1- Computing the roots of S 

We show here that S has all its roots in R[[r/]], by a straightforward application of 
Hensel's lemma. 

First, we need some notation. Given positive integers k,£ and a subset A C R[[??]] fe , 
A mod rf denotes the set [a mod rf \ a £ A]. Besides, we usually denote objects over 
R [[??]] with a * superscript, to distinguish them from their counterparts over R. Finally, 
we extend the notation n to denote the following projection 

tt: R[fo]]» - RN]"" 1 

Recall in what follows that V C R™ is equiprojectable of multi-degree d, that its associ- 
ated triangular set is U, and that S = r/T + (1 — r])U. 

Proposition 5 There exists a unique set V* in R[[?7]] n such that the following holds: 

• V = V* mod 77; 

• V* is equiprojectable of multi-degree d; 

• V* supports interpolation; 

• S is the triangular set associated to V* . 

Proof. We first claim that for i < n and a — (ai, . . . , a n ) in V, the partial derivative 

dS if n , dUi dUi 
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is non zero. Let indeed a' = (ct\, . . . , ct!j_i) £ R* 1 . Then, we have by construction 

Ui(a 1 ,...,a i - 1 ,X i }= Y[(Xi-a), 

so that the previous partial derivative equals 

^^("l) • • • ,ai-i>a») = II («»-«)• 

Since V supports interpolation, this quantity is a product of units, so it is a unit as well, 
establishing our claim. 

Since the system S is triangular, its Jacobian determinant is the product of the partial 
derivatives dSi/dXi. By the previous remark, all these derivatives are units at a, so the 
Jacobian itself is a unit at a. As a consequence, by Hensel's lemma, for all a in V, there 
exists a unique a* in R[[r/]]" such that a = a* mod rj and S(cn*) = 0. We thus let V* C R" 
be the set of all such a* ; hence V = V * mod 77 and S vanishes at V* . 

Next, we prove that V* is equiprojectable of multi-degree d. By induction, we can 
assume that we have proved that tv(V*) is equiprojectable of multi-degree (d\, . . . , d n -\)\ 
it suffices to prove that for each (3* in ir(V*), the fiber Vg* has cardinality d n . 

Let thus a* be in V^„. We prove that for all 7* in V* , ir(a*) = tv(j*) if and only if 
tv (a) — 77(7), with a = a* mod 77 and 7 = 7* mod rj. To prove our claim, remark first 
that if Tv(a*) — tv{^*) then tv{o) = 77(7), by reduction modulo rj. Conversely, suppose 
that tt (a) = 7r(7). Since the system S is triangular, and since a* and 7* are obtained by 
lifting a and 7 using this system, we deduce that n(a*) — 77(7*), as requested. Thus, V* 
is equiprojectable of multi-degree d. 

Finally, we prove that V* supports interpolation. This is again done by induction: 
assume that the projection tv(V*) supports interpolation, let [3* be in tv{V*), and let a* 
and a'* in v^. By assumption on V, a — a' mod i] is a unit in R; thus, by Hensel's lemma, 
a* — a'* is a unit in R[[ry]], as requested. 

This proves the existence of V* with the requested properties. Uniqueness follows in 
a straightforward manner from the uniqueness property of Hensel's lemma. □ 

We continue with complexity estimates: we prove that the roots of S can be computed 
in quasi-linear time. 

Proposition 6 Given T, V and£ > 0, one can compute V* mod rj 1 in time (3(<5 d L(d)M(^)). 

Proof. As before, we proceed inductively: we suppose that the projection W* = tv(V*) 
is known modulo ?/, and show how to deduce V* mod if . To do so, we evaluate all 
coefficients of S n at all points of W* modulo rf . Then, for each (3* in W* , it suffices to 
use Hensel's lemma to lift the roots of S n (/3*,X n ) at precision I. The pseudo-code is in 
Figure 6; for simplicity, we write there W* instead of W* mod rf . 

Lemma 1 shows that we can lift the power series roots of a bivariate polynomial of 
degree d at precision I in time C(d)M(£). As a consequence, the overall cost Cuft of the 
lifting process satisfies 

Cuft(di,.. .,d n ,£) < Cuft(di, . . .,dn-i,£) + C Eva \{di, . . . , d n -i)d n M(£) 
+d 1 ---d n - 1 Q(d n )M(£)- 
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LiftRootsMultivanate(V, S, I) 




1 n=|S| 




2 if n = return [] 




3 W* <- LiftRootsMultivariate(7r(V"), (Si, . . 


,Sn-l),0 


4 for i = 0, . . . , d„ — 1 do 




4.1 valj <- Eval(coeff(S„,X n ,i),W*) 




/* all computations are done modulo rf */ 


5 for /3* in W* do 




5-1 ^-E i<dn val^X n + X^ 




5.2 <- LiftRoots(S /3 *,w /3 ,^) 




6 return [t$„ | /?* e VF*] 





Fig. 6. Lifting the roots of S. 



the middle term gives the cost of evaluating the coefficients of S„ at FF"* mod rf (so we 
apply our evaluation algorithm with power series coefficients); and the right-hand term 
gives the cost of lifting the roots of S„ . This gives 

C Lift (di, . . . ,d n ,£) < ^ C Ev ai(di, • • • , rfi-i)rfiMM + ^ di ■ ■ ■ dj_iC(di)M(*). 
As in the proof of Proposition 3, one deduces that the overall sum is bounded by 

35d c^) M ^ = 3(5dL(d) M ^ ; 

which concludes the proof. □ 

4-2. Multiplication by homotopy 

We continue with the same notation as before. To multiply two polynomials A, B S 
Span(Af d ) modulo (T), we may multiply them modulo (S) over R[rj\ and let r\ = 1 in 
the result. Now the results of the multiplication modulo (S) over R[r/\ and over R[[rj\] are 
the same. When working over R[[f?]], we may use the evaluation-interpolation techniques 
from Subsection 3.5. Indeed, by Proposition 5, S is associated to a subset V* of R[[ry]]™ 
that supports interpolation. 

Of course, when multiplying A and B modulo (S) over R[[?j]], we cannot compute 
with (infinite) power series, but rather with their truncations at a suitable order. On 
the one hand, this order should be larger than the largest degree of a coefficient of the 
multiplication of A and B modulo (S) over R[rj\. On the other hand, this order will 
determine the cost of the multiplication algorithm, so it should be kept to a minimum. 
For e = (ei, . . . , e„) in N™, we define 

H (ei, • • • , e n ) = deg(X 1 ei ■ • • mod (S), v ) 

and 

H(e 1 ,...,e n ) = max ^oK, ■ • • , e'J. 

e^<ei, e' n <e n 

The following proposition relates the cost of our algorithm to the function H; the behavior 
of this function is studied in the next subsection. 
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Proposition 7 Given A, B , T and V , one can compute AB mod (T) in time 0(5dL(d)M( 
0~(5 d r), with r = H(2d 1 - 2, . . . , 2d n - 2). 



Proof. The algorithm is simple: we compute U and use it to obtain V* at a high enough 
precision. In R[Xi, . . . , X n ), the product AB satisfies deg(AB, JQ) < 2di — 2 for all i < n; 
since the multiplication algorithm does not perform any division by n, it suffices to 
apply it with coefficients in R[r]]/ (n r+1 ) , with r = H(2d\ — 2, ... , 2d n — 2). The resulting 
algorithm is given in Figure 7; as before, we write V* for simplicity, whereas we should 
write V* mod rj r+1 . 



Mul(A B, T, V) 

0. U <— AssociatedTriangularSet(V, n) 

1. S <- r/T + (1 - »y)U 

2. <- LiftRootsMultivariate(V r , S, r + 1) 

3. <- Mul(A,B,V*) 

4. return 0,(1, Xi, ... ,X„) 

/* (7,, is seen in R[r]] [Xi , . . . , X n ] */ 



Fig. 7. Multiplication algorithm. 

The computation of U takes time 0(<5dL(d)) by Proposition 3; that of S takes time 
0(<5d). Computing V* mod -q r+1 takes time (3((5dt(d)M(r)) by Proposition 5. Finally, 
the modular multiplication takes time 0((5dF(d)M(r)) by Proposition 4; remark that 
this algorithm is run with coefficients in R[[?7]]/(r7 r+1 ), where all arithmetic operations 
take time <3(M(r)). Finally, specializing r\ at 1 takes time O(Sdr). Summing all these 
costs gives our result. □ 



4-3. Precision analysis 

We show finally how the monomial structure of the polynomials in S affects the cost of 
the algorithm, by means of the integer r of Proposition 7. For i < n and v = (v\, . . . , Vi) 
in N l , we will use the notation = X" 1 ■ ■ ■ Xp and we write the monomial expansion 
of Si as 

s i = xf*+'£ i 8 v x!, (4) 

v£Ei 

where Ei is the set of exponents that appear in Si, the exponents v arc in N 4 , and s v is 
linear in 77. Let us further introduce the coefficients hi defined by ho = and for i > 1, 

h\V\ H h hi-\Vi-i + 1 

hi = max . (5) 

veE t di — Vi 

One easily checks that all hi are positive. The following proposition shows that through 
the coefficients hi, the support Ei determines the cost of our algorithm. 
Proposition 8 The inequality 

H{ex,...,e n ) < h 1 e 1 -\ h 

holds for all (ei, . . . , e n ) £ W l . 
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Using Proposition 7, this proposition gives as an easy corollary the following statement, 
where we take = 2di — 2 < 2di\ we continue using the previous notation T and V. 
Corollary 1 Given A, B, T and V, one can compute AB mod (T) in time 0(5dL(d)M(r)) C 
0~{5 d r), with r < 2(/ii<ii H h h n d n ). 

Hence, the lower the hi the better. However, without putting extra assumptions on the 
monomial supports Ei, Corollary 1 only yields estimates of little interest. Even in sparse 
cases, it remains difficult to simplify the recurrence giving the coefficients hi. Still, several 
examples in the next section will show that for some useful families of monomial supports, 
significantly sharper bounds can be derived. 

The rest of this section is devoted to prove Proposition 8. In all that follows, the 
multi-degree d = (di, . . . , d n ) and the supports Ei are fixed. We also let E[ be the set of 
modified exponents 

E\ = {v - (0, . . . , 0, di) e Z* | v e Ei}, 
so that for all v = (ui, . . . , Vi) in E\, Vj > for j < i and z/, < 0. Hence, Equation (5) 
takes the (slightly more handy) form 

h\V\ H V hi-m-! + 1 

hi = max . (6) 

v£E'. —Vi 

Recall that the function Hq was defined in the previous subsection with domain N"; 
in what follows, we also consider Ho as a function over W, for 1 < i < n, by defining 
Ho(ei, . . . , ei) = -ffo(ei, . . . , ej, 0, . . . , 0), where the right-hand expression contains n — i 
zeros; for completeness, we write Hq() = for i — 0. The following recurrence relation 
enables us to control the growth of H . 

Lemma 2 For i > 1, let e = (ei, . . .,ej) be in W and let e' — (ei, . . . ,ej_i) in N 1 ^ 1 . 
Then the following (in) equalities hold: 

Ho(e) = Ho(e) if e-i < di, Ho(e) < 1 + max Ho(e + v) otherwise. 



Proof. Let us first suppose a < di] then, 

■ ■ ■ X? mod (S) = • • • mod (S»X?', 

since the latter product is reduced modulo (S). Both sides have thus the same degree in 
77, and our first claim follows. 

We can now focus on the case > d,, for which we write /j = — di, so that fi > 0. 
From Equation (4), we deduce 

X{*Si = x{* +d > + a^/'Xr, 

and thus we get 

X f i >S i =X?+Y J s v X^\, 
veE> 

by the definition of E[. In our notation, we have X^ 1 • • • Xf t l — Xf . Thus, after multipli- 
cation by X^ 1 • • • X^ 1 and term reorganization, the former equality implies that 

x? - • • • xf^x^Si = - s ^t +v - 
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As a consequence, we deduce that 

deg(X* mod (S),r?) < maxdeg^Xf^ mod (S),?7). 

v£E[ 

Since for v in E[, we have 

deg( S „X^ mod (S), n) = degK, r?) + deg(Xf +y mod (S), r?) = 1 + H (e + v), 
the conclusion follows. □ 

Iterating the process of the previous lemma, we obtain the following bound. In the 
next lemma, (f u )veE' are a family of integer valued variables. 
Lemma 3 Let e = (ei, . . . , ei) be in W . Then the following inequality holds: 

flo(e) < max H (e + £ f v v) + £ /„. 

{fv)v£E'. non-negative integers "£ E l V ^ E [ 

such that < ei + 5Z„ gs ; f v Vi < d» — 1 

Proof. We prove the claim by induction on e^. For ei < di — 1, the family (/„ = 0)„ £ £' 
satisfies the constraint < ei + Y^veE' fv v i < ^ — 1; for this choice, the value of the 
function we maximize is precisely H (e), so our claim holds. Suppose now that ei > di. 
Then, the previous lemma gives 

H (e) < 1 + max H (e + v). (7) 

v&E'. 

Let us fix v va E'i, then e + z/ has non-negative integer coordinates, and its ith coordinate 
is less than a. Thus, we can apply the induction assumption, obtaining 

H (e+i>) < max H {e+v+ ^ f„>v')+ E U> ■ 

{fv')v'eE'. non-negative integers v ' eE 'i u ' eE 'i 

such that < a + ut + J2„>£E' f^' v i < di — I 

To any set of non-negative integers (fu')v'eE' with 

0<ei + Ui+ !w[<d l -\ 
v'eE 1 . 

appearing in the previous maximum, we associate the non-negative integers (fli)u'eE', 
with f v = f v + 1 and f v , = f v i otherwise. These new integers satisfy 

0<e t + f>'i<di-l 
v'eE- 

and 

H a {e + u+ fv'"')+ E f'=H (e+ E />') + E ~ L 

v'eE[ V '£- E 'i v '£ E 'i v '^ E 'i 
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Taking maxima, we deduce from the previous inequality 

H (e + u)< max H Q (e + ^ />') + £ f v , - 1. 

(fl')v'eE'. non-negative integers i>'e.EJ 
such that < e, + ^„, eB / < di — 1 

Substituting in Equation (7) and taking the maximum over v in E\ concludes the 
proof. □ 

For i < n, let Li be the linear form (ei, . . . , e*) i— > /iiei + • • • + /i^i, where the ftj are 
as in Equation (5). The following lemma concludes the proof of Proposition 8; as we did 
for H , for i < n, we extend H to N l , by writing H(e±, . . . , ej) = H(ei, . . . , ej, 0, . . . , 0). 
Lemma 4 For i < n and e = (ei, . . . , a) in W , the inequality H(e) < Lj(e) holds. 

Proof. It is sufficient to prove that H (e) < Lj(e) holds; since all coefficients of Li are 
non-negative, Li is non-decreasing with respect to all of its variables, which implies the 
thesis. 

We prove our inequalities by induction on i > 0. For i = 0, we have -ffoO = ^o() = 0; 
hence, our claim vacuously holds at this index. For i > 1, we now prove that if our 
inequality holds at index i — 1, it will also hold at index i. Lemma 3 shows that for any 
eef , we have the inequality 

H (e) < max H (e+ ^ f v v) + £ /„. 

(/„)„ eE < non-negative integers U ^ E [ V ^ E ', 

such that < e» + E„ 6 £' /f < d 4 — 1 

Let V? be the natural projection N 1 — ► N' _1 , let {j v ) ve E' be non-negative integers that 
satisfy the conditions in the previous inequality. Since e + 5Z„ e£ v /i/f has degree in 
less than d i; the first point of Lemma 2 shows that 

the induction assumption implies that this quantity is bounded from above by 

L i _ 1 (ip(e+ E ./»)• 

As a consequence, H (e) admits the upper bound 

max L i - 1 (tf(e+ ^ U v )) + E -f"- 

{fv)veE'. non-negative integers V ^ E [ V ^ E [ 

such that < ei + $^„ gB ; f v Vi < di ~ 1 

i 

This quantity itself is upper-bounded by a similar expression, where we allow the f v to 
be non-negative reals numbers; this gives 

H (e) < max L i _ 1 ((p(e + ^ f» v ))+ E 

(/„)„ eE / non-negative real numbers v ^ E 'i V ^ E [ 

such that < e, + X^es' /i/^; < dj — 1 
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Since all hi and all v\, . . . , Vi-\ arc non-negative, the function of {f u )veE' wc want to 
maximize is affine with non-negative coefficients. The domain where we maximize it is 
defined by the conditions 

/„ > for all v € E[, < a + ^ f v v { < d { - 1, 

and it is contained in the domain D defined by the conditions 

f v > for all v e El, < e % + f" u i- 

v eE> 

Since all unknowns /„ are non-negative, while the coefficients Vi are negative, the domain 
D is convex and bounded. Hence, the maximal value we look for is upper-bounded by 
the maximal value at the end-vertices of D, distinct from the origin; these vertices are 

E v = {f v ,= for v' ± v, = for v e E\. 

Vi 

At the point E v , the objective function takes the value 

Li_i(^(e- — v)) - —. 

Vi Vi 

By the linearity of Li_i and ip, this can be rewritten as 

Li_i(<p(e)) - Li_i(ip(— v)) = Li_i((p(e)) e { . 

Vi V^ Vi 

As a consequence, we obtain the upper bound 

H (e) <£ t _ 1 ( y (e)) + max Z/ - l(y(t/)) + 1 et . 

veE\ -Vi 

To simplify this further, note that the term L i _ 1 ((/3(e)) rewrites as h Y e\ + • • • + /ij_iej_i. 
Similarly, Li-i(ip{v)) + 1 equals h\V\ + ■ ■ ■ + hi-\Vi_\. We deduce the inequality 

rr / w , . , ^\V\ -\ h hi-XVi-X + 1 

H (e) < h\e\ + ■ ■ ■ + hi-\ei-\ + max a, 

>'C/;; -//; 

which we can finally rewrite as 

H (e) < h\e\ H h /ij_iej_i + ft^i, 

as requested. □ 



5. Examples 

5. 1 . Main family of examples 

We give explicit estimates for the coefficients hi of the previous section on the following 
family of examples. We consider triangular sets T = (T\, . . . ,T„) such that Tj has the 
form 

Ti = X* + ]T t v X?, t„ G R, (8) 
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where all are monomials in X\, . . . ,Xi of total degree at most Aj, for some Aj € N. 
We let d — maxi<„ di, and we suppose that R contains at least d pairwise distinct values 
Xi,. . . ,x d , with Xi — Xj a unit for i ^ j. 

The following proposition illustrates three different situations. The first two cases 
display a cost quasi- linear in dd d , which is satisfying, especially for small d; the last one 
shows that small changes in the assumptions can induce large overheads. We will sec in 
the next subsection cases where di is constant equal to d, or di = n+ 1 — i; in such cases, 
d is logarithmic in 5 d and the cost 0~(d5 d ) is thus 0~(5 d ), which is what we were aiming 
at. 

Proposition 9 With assumptions as above, multiplication modulo (T) can be performed 
with the following complexities: 

0(nS d ^M(nd)) c 0~(dS d ) if Xi = di - 1 for all i, 

0(n5 d ^f-U{n 2 d)) c 0~(dS d ) if A* = di for all i, 

0(nS d ^p-M{2 n d)) C 0~(2 n dS d ) if X t = d { + 1 for all i. 



Proof. First, we construct V: we simply choose the grid 

V = x ••• x [x\, . . . ,x dn ]. (9) 

Thus, we have Ui = (X t - zi) • • • (X, - x di ); as before we let S = rfT + (1 - »j)U. Thus, 
the monomial support Ei associated with Si is contained in 

D' i = £>jU{(0,...,0,i/i) \Q<Vi<di}. 

Since each monomial in Di has an exponent of the form (y\, . . . , fj), with V\ + • • • + Vi < Aj 
and Vi < di, we deduce from Equation (5) that 

ftii/H hfti-i^i-i + l / ftii/H h /ij-ifj-i + 1 A 

/ij < max < max max , 1 1 . 

di — Vi V ^efli di — Vi J 

Let /i^ = max(/ii, . . . , /ij), so that 

. . / + + + l A . / &i_i(Aj-l/j) + l \ 
«i < max max ,1 < max max , 1 1 . 

\veDi di — Vi ) \veDi di — Vi ) 

(10) 

Knowing the distribution of the di and A^, the former relation makes it possible to analyze 
the growth of the coefficients hi, and thus of 2(d 1 /i 1 + • • • + d n h n ). 

Case 1. Suppose first that Xi = di — 1. Then, the former inequality implies hi < 1 for 
all i, so that 2{d\hx + • • • + d n h n ) < 2nd. 

Case 2. If Aj = di, then (10) becomes hi < h' i _ 1 + 1, so that hi < i for all i, and thus 
2{d\h\ + • • • + d n h n ) < n{n + l)d. 

Case 3. If finally Xi = di + 1, then (10) becomes hi < 2h' i _ 1 + 1, so that h\ < 2* — 1. 
In this case, we get 2{d\hx + • • • + d n h n ) < 2 n+2 d. 

To conclude the proof, we simply plug the previous estimates in the cost estimate 
0(S d L(d)M(r)) of Corollary 1, with r < 2(d\hi + • • • + d n h n ), and we use the upper 
bound L(d) < nC{d)/d of Equation (3). □ 
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5.2. Cauchy modules 



Cauchy modules [37] are a basic construction in Galois theory and invariant theory [32 
[271[TJ [55]. Starting from a monic polynomial F £ R[X] of degree d, we define a triangular 
set Fi,...,Fd by letting F 1 (X 1 ) = F(Xi) and taking iterated divided differences: 

c, fY Y n Fi(Xu ■ ■ ■ >Xi-i,Xi) - Fi(X lt . . . ,Xi-i,X i+1 ) 

fii+x{Xi, . . . ,x i+ i) = — — i < i < a. 

The polynomials Fi, . . . , Fd form a triangular set of multi-degree d = (d, d— 1, . . . , 1), so 
that 5d — dl; their interest stems from the fact that they form a system of generators of 
the ideal (cr, — (— l) 1 fd-i)i<i<d, where o~i is the ith elementary symmetric polynomial in 
Xi, . . . , Xd and fi is the coefficient of X 1 in F. 

One easily checks that Fi has total degree at most d + 1 — i. Hence, assuming that 
0, . . . , d— 1 are units in R, we are under the assumptions of Subsection 5.1, with A; = di = 
d + 1 — i for alii and (xi, . . . , Xd) = (0, . . . , d— 1). As a consequence, Proposition 9 shows 
that multiplication modulo (F±, . . . ,Fd) can be done using 0(d\C(d) M(d 3 )) operations 
in R, that is, in quasi-linear time 0~(dl). This improves for instance the results given 
in |17j on the evaluation properties of symmetric polynomials. 

5.3. Polynomial multiplication 

We show now how to derive quasi-linear time algorithms for univariate multiplication 
in R[X] from our previous multivariate construction. Unfortunately, our algorithm does 
not improve on the complexity of Cantor-Kaltofen's algorithm 9 ; however, we believe 
it is worth mentioning. Precisely, given n > 1, we give here an algorithm to perform 
truncated multiplication in R[X]/(X 2 ). We introduce variables X\, . . . ,X n ; computing 
in A = R[X]/ (X 2 ) is equivalent to computing in B = R[Xi, . . . , X n )/(Vi, . . . , V n ), with 
V= (V u ..., V n ) given by 

Xi - X n 
X n -i — X 2 

since the dummy variables X\, . . . , JC n _i play no role in this representation. However, 
changing the order of the variables, we see that the ideal (Vx, . . . , V n ) is also equal to the 
ideal (Ti, . . . ,T n ) given by 

X 2 — X n -i 
X 2 - Xi 

xi 

The R-basis of B corresponding to V is (X^)i<2»; the basis corresponding to T is Md 
(notation defined in the introduction), with d = (2, . . . , 2). Besides, the change of basis 
does not use any arithmetic operation, since it amounts to rewrite the exponents i in 
base 2, and conversely. 
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Hence, we can apply our multivariate multiplication algorithm modulo (T). Remark 
that the triangular set T satisfies the assumptions of Subsection 5.1 (for any R), with d\ = 
■ ■ ■ = d n = d = 2, Sd = 2 n , Ai = • • • = A„ = 1 and (0:1,2:2) = (0, 1). By Proposition 9, 
we deduce that the cost of a multiplication in B, and thus in A, is 0(2 n nM(n)). Since 
one can multiply univariate polynomials of degree 2™ using two multiplications in A, this 
gives the recurrence 

M(2 n ) < k2"nM(n) and thus M(d) < k'tf log(d)M(log(d)) 

for some constants k, V! . Unrolling the recursion 1,2,..., times, and taking M(n) £ 0(n 2 ) 
to end the recursion, we obtain quasi-linear estimates of the form 

M(d) e 0(dlog(d) 3 ) or M(d) E 0(dlog(d) 2 log(log(d)) 3 ), . . . 

The main noteworthy feature of this multiplication algorithm is that no root of unity is 
present, though our multivariate evaluation-interpolation routine is somewhat similar to 
a multivariate Fourier Transform. In particular, the case when 2 is a zero-divisor in R 
requires no special treatment, contrary to [S]- 

5.4- Exponential generating series multiplication 

We continue with a question somehow similar to the one in the previous subsection. 
Given two sequences do, • • • , a>d and bo, ■ ■ ■ ,bd in R, we want to compute the sequence 
Cq , . . . , Q such that 




where the binomial coefficients are the coefficients of the expansion of (1 + X) 1 in R[X]. 
We discuss an application of this question in the next section. 

The naive algorithm has cost 0(d 2 ). If 1, ...,d are units in R, the former equation 
takes the form 

i<d i<d i<d 

so we can achieve a cost 0(M(d)). Under some much milder assumptions on R, we are 
going to see how to achieve a similar cost through multivariate computations. 

We will suppose that there exists a prime p such that for a £ N, if gcd(a,p) = 1, then 
a is a unit in R (this is the case e.g. for R = Z/p k Z). Let n be such that d+1 < p n , and 
introduce the triangular set T = (Ti, . . . , T n ) defined by 

XP-pX^ 

xl- v x x 
xf. 

In what follows, for i > 0, (io, i\, . . . ) denotes the sequence of its coefficients in base p; 
thus, for i < d, only %q, . . , , i n -i can be non-zero. Besides, we let / : N — > N be defined 
by fii) = i\/p vtyl '\ where v(il) is the p-adic valuation of i\. In particular, f(i) is a unit 
in R. 



21 



Proposition 10 Let 



i<d J w i<d J v ' i<d J y ' 

Then C = AB mod (T) . 

Proof. Let i,j<d, with k = i + j < d. We start by the obvious remark that 

k\ _/(*)_»((*)) (13) 

holds in R. Besides, the normal form of the product -jj^X l ° ■ ■ ■ Aj" -1 by j^X^ ■ ■ ■ X^"^ 1 
modulo (T) is 

a ' b J n c(i,j) x k ° ■■■ X kn - X 

mm* - 1 ' 

where Ci.j is the number of carries held in the addition of i and j in base p. From [31 L 
Eq. (1.6)], Cij is exactly the valuation of the binomial coefficient Thus, by (13), the 
former product equals 



l 



,j xl° ■■■ x\' 



Summing over all i,j gives our claim. □ 

As in the previous subsection, we can apply our multivariate multiplication algorithm 
modulo (T). Remark that the triangular set T satisfies the assumptions of Subsection 5.1, 
with di = ■ ■ ■ = d n = d = p, 6^ = p n , Ai = • • • = A n = 1 and (x%, . . . , x p ) = (0, . . . ,p— 1). 
Note as well that we can take n S 0(log p (d)), and that 5^ — p n < pd. 

By Proposition 9, we deduce that the cost of computing C, and thus all Cq, . . . , c<j, is 
0(dlog(d) M(p) M(plog p (d))). If p is fixed, we obtain the estimate 0(dlog(d)M(log(d))). 
This is not as good as the estimate 0(M(d)) = 0(dlog(d) loglog(rf)) we obtained in 
characteristic zero, but quite close. 



6. Application: computing with algebraic numbers 

We finally present an application of the previous constructions to computation with 
algebraic numbers, and give timings of our implementation. 

6.1. Presentation of the problem 

Let k be a field and let / and g be monic polynomials in k[T), of degrees m and n 
respectively. We are interested in computing their composed sum h = / ® g. This is the 
polynomial of degree d — mn defined by 

/® 5 = IJ(T-a-/3), 

the product running over all the roots a of / and /3 of g, counted with multiplicities, in 
an algebraic closure k of k. 
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A natural approach consists in computing h(T) as the resultant of f(T — U) and g(U) 
in U. However, the fastest algorithm for resultants [25 j has a complexity of order 0~(d 15 ) 
for m = n. To do better, Dvornicich and Traverso [12] suggested to compute the power 
sums 

f(a)=0 9(/9)=0 

of respectively / and g, and deduce the power sums Cj of h, by means of Equation (11). 
In [5J, this approach is showed to take time 0(M(d)), over fields of characteristic zero 
or larger than d. Indeed, computing (di)i<d and (6j)j<d can be done in 0(M(d)) oper- 
ations, over any field, using Newton iteration for power series division [28 . Then, by 
our assumption on the characteristic, one can compute (cj)j<d in quasi-linear time using 
Equation (12), for another M(d) + O(d) operations. Finally, knowing (ci)i<d, one can 
then recover h in time 0(M(d)) as well, using fast exponential computation [51 12"51 IMl IT] ; 
this step relies as well on the assumption on the characteristic. 

If k has positive characteristic less than d, two issues arise: Equation (12) makes no 
sense anymore and (ci)i<d are actually not enough to recover h. To our knowledge, no 
general solution better than the resultant method was known up to now (partial answers 
are in [51 129j under restrictive conditions). We propose here a solution that works over 
finite fields, following an idea introduced in [T5] , 

For simplicity, we consider only k — ¥ p . Since our algorithm actually does computa- 
tions over rings of the form Z/p"Z, measuring its complexity in F p -operations as we did 
up to now is not appropriate. Instead, we count bit operations. Thus, we let Mz be such 
that integers of bit-length t can be multiplied using M z (£) bit operations; quasi-linear 
estimates are known as well for M z , the best to date being Fiirer's £log(£)2°( log ^ [14] . 
Proposition 11 Given f and g, one can compute h using 

0((M(d) + d log(d) M(p) M(plog d (p))) N(p,d)) 

bit operations, with N(p,d) = 0(M z (log(p)) log(log(p)) + M z (log(d))). 

After simplification, this cost is seen to be 0~(dp 2 ) bit operations. Also, if we consider p 

fixed, the cost becomes 

0((M(d) + dlog(d)M(log(d)))M z (log(d))), 

that is, quasi-linear. 

Proof. Let Z p be the ring of p-adic integers and let F and G be monic lifts of / and 
g in Z p [T], of degrees m and n. Defining H = F G £ Z p [T], we have that h — 
H mod p. Let further (Ai)i>o, (-Bi)i>o and (Ci)i>o be the power sums of respectively F, 
G and H. For any a > 0, the reductions Ai mod p a , Bi mod p a , and Ci mod p a satisfy 
Equation (11), so we can apply the results of Subsection 5.4 to deduce (Ci mod p a )i<d 
from (A4 mod p a )i< d and (B t mod p a )i< d - 

Besides, taking a = |log p (d)J + 1, it is proved in jS] that given (Ci mod p a )i<d, one 
can compute h in quasi-linear time 0(M(d)M z (log (d))) bit operations. Remark that this 
step is non trivial: recovering a polynomial of degree d from its Newton sums requires 
divisions by 1, . . . , d, and not all these numbers are units in small characteristic. 

In the algorithm, the function Lift simply lifts its argument from ¥ P [T] = Z/pZ[T] to 
Z/p a Z[T}; the function PowerSums computes the first d power sums of its arguments by 
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the algorithm of [55]. Step 7 applies the algorithm of Subsection 5.4, and the last step 
uses the algorithm presented in [B] to recover h. 

Our choice of a implies that log(p a ) = 0(log(rf)). Thus, operations (+, x) modulo p a 
take 0(Mz(log(rf))) bit operations [15| Chapter 9]. Using Newton iteration, inversions 
modulo p a take N(p, d) = 0(Mz(log(p)) log(log(p)) + Mz(log(d))) bit operations, where 
the first term stands for the cost computing the inverse modulo p, and the second one 
for lifting it modulo p a . 

The cost of computing {Ai)i<d and (i?i)i<<j is 0(M(d)) operations modulo p a ; this 
dominates the cost of recovering h. The remaining cost is that of computing (Ci)i<d, 
which is reported in Subsection 5.4 in terms of numbers of operations modulo p a . The 
previous estimate on N (p, d) concludes the proof. □ 



ComposedSum(/, g) 

1. d <- deg(/) deg(g) 

2. a «- Llog p (d)J + 1 

3. F <— Lift(/, a) 

4. {Ai)i<d <— PowerSums(.F, d) 

5. G^~L\ft(g,a) 

6. (Bi)i<d <— PowerSums(G, d) 

7- (Ci)i<d *— ExponentialGeneratingSeriesMultiplication(A, i?) 

8. return PowerSumsToPolynomial(C) 



Fig. 8. Composed sum in small characteristic. 



6.2. Experimental results 

We implemented the composed sum algorithm over F 2 (i.e., p — 2 here). We used the 
NTL C++ package as a basis [30]. Since NTL does not implement bivariate resultants, we 
also used Magma [3] for comparison with the resultant method. All timings are obtained 
on an AMD Athlon 64 with 5GB of RAM. 

Figure 9 gives detailed timings for our algorithm; each colored area gives the time of 
one of the main tasks. The less costly step is the first, the conversion from the original 
polynomials to their Newton sums. Then, we give the time needed to compute all the 
power series roots needed for our multiplication algorithm, followed by the evaluation- 
interpolation process itself; finally, we give the time necessary to recover h from its 
power sums. Altogether, the practical behavior of our algorithm matches the quasi-linear 
complexity estimates. The steps we observe correspond to the increase in the number of 
variables in our multivariate polynomials, and are the analogues of the steps observed in 
classical FFT. 

Figure 10 gives timings obtained in Magma, using the built-in resultant function, on 
the same set of problems as above. As predicted by the complexity analysis, the results 
are significantly slower (about two orders of magnitude for the larger problems). 
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7. Conclusion 



Several questions remain open after this work. Of course, the most challenging one 
remains how to unconditionally get rid of all exponential factors in multiplication algo- 
rithms for triangular sets. More immediate questions may be the following: at the fine 
tuning level, adapting the idea of the Truncated Fourier Transform [33] should enable 
us to reduce the step effect in the timings of the previous section. Besides, it will be 
worthwhile to investigate what other applications can be dealt with using the "homo- 
topy multiplication" model, such as the product of matrices with entries defined modulo 
a triangular set, or further tasks such as modular inversion or modular composition. 
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Fig. 9. Detailed timings for our algorithm 
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Fig. 10. Timings in magma 
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